Search for non-Gaussianity in pixel, harmonic and wavelet space: 

compared and combined 

Paolo Cabella/'Q Frode Hansen/ Domenico 



Marinucci 



; 2 







Via della Ricerca Scientifica 1, 1-00133 Roma, Raly 



Daniele Pagano/ and Nicola Vittorio^ 

^ Dipartimento di Fisica, Universitd di Roma 'Tor Vergata', 

Via della Ricerca Scientifica 1, 1-00133 Roma, Raly 

^ I ^ Dipartimento di Matematica, Universitd di Roma 'Tor Vergata', 

(N 

^ ! (Dated: February 2, 2008) 

^ Abstract 

We present a comparison between three approaches to test non-Gaussianity of cosmic microwave 
' background data. The Minkowski functionals, the empirical process method and the skewness of 

o ■ 

' wavelet coefficients are applied to maps generated from non-standard inflationary models and to 

. Gaussian maps with point sources included. We discuss the different power of the pixel, harmonic 

'NT 

o 

Qh| also suggest a new procedure consisting of a combination of statistics in pixel, harmonic and wavelet 

o 

■ space. 

to 



and wavelet space methods on these simulated almost full-sky data (with Planck like noise). We 



PACS numbers: 02.50.Ng, 95.75.Pq, 02.50.Tt, 98.80.Es 



* Electronic address: 'Paolo. Uabella(y)rom a2.infn.itl 
Electronic address: frodekh@roma2.inf n.'itl 
•l-Electronic address: lmarinucc@mat . uniroma2ltl 



1 



I. INTRODUCTION 



The fluctuations of the cosmic microwave background (CMB) are expected to be close 
to Gaussian distributed. In view of the increasing quantity of CMB experiments, it is 
now possible to check this assumption on data with growing resolution and sky coverage. 
Most models for the early universe predict some small deviations from Gaussianity; 



ersepreaict some s: 



n 



non-standard models of inflation [1|, 12, |3|, y, |5|, |6|, |7[, cosmic strings (See Ref. [8| for a 

review) and point sources. Detecting these small deviations would be of great importance 
for the understanding of the physics of the early universe. Also systematic effects like a 
non-symmetric beam and noise could give rise to non-Gaussian features. For this reason a 
non- Gaussianity check could reveal whether the impact of the instrumental effects on the 
data of the experiment is well understood. 



The methods to search for non-Gaussianity in the literature mainly concentrate on im- 

Elementing the test in three different spaces: (1) In pixel space: the Minkowski functionals 
], (which were used to set limits on the non-Gaussianity in the WMAP data 
temperature correlation functions the peak to peak correlation function 1^, skewness 
and kurtosis of the temperature field [l^ and curvature properties jl^ Il6|, to mention 
a few. (2) In harmonic space: analysis of the bispectrum and its normalized version 
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201 1 and the bispectrum in the flat sky approximation j21[. The explicit form of 



the trispectrum for CMB data was derived in 



22l . l23l | . Phase mapping |2J| . Applications 



to COBE, Maxima and Boomerang data have also drawn enormous attention and raised 



wide de 



debate J25j26|j27 
: 3iH34,|3j. 
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3l|. Finally, (3) wavelet 



28j . The empirical process method 
space: |^ yj, |^J. Traditionally, these tests are performed separately in each space. 
In this article, we will take methods in pixel- (the Minkowski functionals), harmonic- (the 
empirical process) and wavelet-space (skewness), and we will make a comparison for two 
different models of non-Gaussianity. We will also combine the methods in order to improve 
the total power. It should be noted that all the procedures we consider are non-parametric, 
that is they do not assume any a priori knowledge about the nature of non-Gaussianity. 



We will use these methods on 100 maps generated from a non-standard inflationary 
model 0] and on 100 maps where we have included point sources. We assess the performance 
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of the methods in the different spaces for the different types of non-Gaussianity. We also 
propose a combined test which turns out to be more robust. 

In sectionini we review the method of Minkowski functionals, in section lTTTl we describe our 
implementation of the method and in section HVI we define our proposed statistic. In section 
IVl we review the empirical process method while section IVTl is devoted to the wavelets. In 
section IVIII the methods are compared and applied to non-standard inflationary models, in 
section IVIIII to maps with point sources. Finally in section HXl we summarize and comment 
on our results. 



II. MINKOWSKI FUNCTIONALS 



To analyze a spherical map in terms of Minkowski functionals, we consider the excursion 
sets, that is, the map subsets which exceed a given threshold value. The threshold is 
labelled u, and it is treated as an independent variable, on which these functionals depend. 
More precisely, considering the normalized random field of temperature fluctuations, u = 
AT/a(AT); we can define the 'hot region' Q as the ensemble of pixels Ui higher than the u 
level : 

Q = Qiu) = {t\u{9„ip,)>u} . (1) 

The three functionals of interest then are, up to constant factors 0|: 

1) Area: Mo(z/) is the total area of all hot regions. 

2) Boundary length: Mi(z/) is proportional to the total length of the boundary between 
cold and hot regions 

3) Euler characteristic or genus: M2(z/), a purely topological quantity, counts the num- 
ber of isolated hot regions minus the number of isolated cold regions, i.e. the number of 
connected components in Q minus the number of 'holes'. 

The rationale behind these statistics can be explained from mathematical results in Had- 
wiger (1959); in particular, these results can be interpreted by stating that all the morpho- 
logical information of a convex body is contained in the Minkowski functionals (Winitzki 
and Kosowsky, 1997)); here, by morphological we mean the properties which are invariant 
under translations and rotations and which are additive The three statistics, 
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normalized by the area density, can then be expressed as 

MM = J [ 

J Q 

M2M = I Kdl 



2nA 



(2) 
(3) 
(4) 



dQ 



where dQ is the contour of the region Q; dA and dl are the differential elements of Q and 
dQ, respectively; k is the geodetic curvature of dl . 

The expected values for a given thresholds depends on a single parameter r given for a 
Gaussian field by 39^ : 



Mo(z/) 

MM 
M2{v) 



1 
2 



1 — erf 



exp 



-.V exp 



with 



1 




r 



(5) 

(6) 
(7) 

(8) 



where semicolon indicates the covariant derivative on the sphere. 
In the case of CMB (IHl) reduces to 



T = ^(2/ + 



/(/ + !) 



(9) 



where Ci is the angular power spectrum. 

An immediate consequence of the above formulae is that, although the expected value 
of the first Minkowski functional is invariant with respect to the dependence structure of 
AT, for the second and third Minkowski functional this is not the case and calibration 
for a given angular power spectrum Ce is needed. Moreover, even for the first Minkowski 
functional, knowledge of the angular power spectrum is required for a Monte Carlo evaluation 
of its variance. This can be viewed as a drawback, and because of this some effort has 
been undertaken to provide at least some crude upper bound for the functionals' variance 
(Winitzki and Kosowsky (1997)). 
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III. THE IMPLEMENTATION OF MINKOWSKI FUNCTIONALS 



In order to estimate the three functionals we simulate a map of the CMB j40[ with a known 
power spectrum, and then we cut the maximum number of independent tangent planes of 
dimension ~ 12° x 12°; in this way it is easy to calculate the values of the three Minkowski 
functionals in the flat-sky limit, taking into account the possibility of gaps (galactic cut, 
polar calottes). Due to projection effects, finite pixel size and the dimension of the tangent 
planes, we find a deviation of the simulated values with respect to the analytical spherical 
predictions (eqs|Hl). Note however that the shape of the curves is unaffected, which is not 
the case when non-Gaussianities are present (see fig. |2I). In fig. Qwe show a comparison 
between analytical expectation values and the values computed on the tangent planes. 

IV. TEST OF NON-GAUSSIANITY IN PIXEL SPACE 

In order to test non-Gaussianity we use a test defined by: 

\M\v) - M\v)\dv i = 0, 1, 2 (10) 



Our procedure is as follows; 

• Given an observed map, we estimate the power spectrum. 

• Using the estimated power spectrum, we generate 200 Gaussian realizations 







For each map we calculate the Minkowski functionals, using the tangent planes as 
described above 

We calibrate the quantiles (i.e. the threshold values at a given significance level) using 
the Monte Carlo simulations; the one and two sigma detection levels are shown in 
figO 

We calculate our statistic (fTUI) from the observed map and compare the result with 
the Monte-Carlo calibration; we are thus able to determine at which confidence level 
the map is Gaussian or not. 
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FIG. 1: Comparison between Minkowski functionals computed in the tangent plane approximation 
and their analytical predictions (solid line) see Eqs.©. Squares refer to tangent planes of size 
~ (12*^ X 12*^), triangles refer to tangent planes of size ~ (24*^ x 24'^). 
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V. THE EMPIRICAL PROCESS METHOD 



The details of the empirical process approach to detect non-Gaussianity in the CMB were 
Efiven m 0,0,0. In short, the method consists of a family of tests which focus on the 
total distribution of aem and check for dependencies between k ^-rows. The first step is to 
transform the spherical harmonic coefficients into variables uem which have an approximate 
uniform distribution between and 1, given that the a£m were initially Gaussian distributed. 
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FIG. 2: A comparison between a Gaussian map (solid line) and non-Gaussian map with fNL = 1000 
(squares) (the /nl factor will be explained in section IVTH) . 

This is done using the Smirnov transformation, defined as 

(I (^£0 P \ Z' 2 1 Ct£m I ^ \ 

^-^j ,uem = ^2\^ ^ j , m = 1,2,...,/, / = 1,2,...L , 

where is the cumulative distribution function of a with n degrees of freedom and 
Ci are the power spectrum coefficients estimated from the data. The error introduced by 
using estimated (7^ instead of the real underlying Ci is dealt with using a bias-subtraction, 
as described in 12911. 
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FIG. 3: Histogram with threshold levels for all three Minkowski functionals for Gaussian realiza- 
tions (solid line) and non-Gaussian realizations with /nl = 100 (dotted line) (the /nl factor will 
be explained in section Wll'^ . All realizations have Planck-like noise and beam 20'. The shaded 
areas represent the la and 2a detection limits. 

Then the joint empirical distribution function for row i is formed, 

where j determines the spacing between the rows for which the dependencies are tested 
and Ami denotes the difference in m for row i. The parameters run over the interval [0, 1]. 
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FIG. 4: Minkowski functionals averaged over 200 realization of Gaussian maps (solid line) and 100 
non-Gaussian maps (squares) with f^L = 300, Planck noise and beam 20'. The small la deviations 
are also shown. 

The empirical process is expressed using the centered and rescaled F^^^i+Ai ^-i given as 

I i=l 

The intuition behind this procedure is as follows: if the a^^s are Gaussian, G converges to 
a well-defined limiting process, whose distribution can be readily tabulated. On the other 
hand, for non-Gaussian a^^s {F{ai, ak) — Y[i=i ^i} thereby G will take 'high' values 
over some parts of a-space. Thus, the analysis of some appropriate functional of G can be 
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FIG. 5: Minkowski functionals (divided by the mean) averaged over 200 realization of Gaussian 

maps (solid line) and 100 non-Gaussian maps (squares) with /jv^ = 300/100, Planck noise and 

beam 20'. The shaded bands show the la deviations of the Gaussian realizations and the error 

bars show the la deviations for the non-Gaussian maps. 
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used to detect non-Gaussianity. To combine the information over all multipoles into one 
statistic, we define 

1 [(i-Afc_i)r] ^ 

KL{ai,...,ak,r) = ^= V G^,...,^+a, ^.^(ai, Ofc) , (11) 

a/-L - iXe^k-i 

where L is the highest multipole where the data is signal dominated. 

The method can then be summarized as follows: the distribution of supli^^^l is found 
using Monte-Carlo simulations of Gaussian distributed aum- Then, for a given observed 
set of airni the value fcmax = suplf^il is found and compared to the distribution obtained 
from Monte-Carlo. The consistency of the data with a Gaussian distribution can then be 
estimated to any suitable a-level. In this simple approach was extended in three different 
ways. First of all, the fact that the above explained estimator is not rotationally invariant is 
exploited, using the fcmax value averaged over many rotations. Each rotation can be viewed 
as a resampling of the a^rn- Secondly, we introduced three variations of the test, taking 
into account, not only the modulus of the a^m but also the phases. Finally, experimental 
effects like noise and galactic cut was accounted for using Monte Carlo calibration of the 
^max distribution with these effects included. It should be noted that the rotated maps are 
clearly dependent and the resulting statistic may thus depend on the shape of the angular 
power spectrum. 



VI. TEST OF NON-GAUSSIANITY IN WAVELET SPACE 



A third space where one could look for non-Gaussianity is the wavelet space. The use 



of wave^ 

y, 33, 



ets for non-Gaussianity tests of the CMB has been investigated by several authors 



35| and turns out to be a very powerful tool. We will here just briefly describe 



the wavelet method, and refer to the above references for more details. 
An isotropic wavelet can be deflned as 



with the properties 



/ 



dx ip{x) 



< oo 



(12) 

(13) 
(14) 
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where x = R represents a scale and b a translation. The Fourier transform of the 
wavelet is represented by ipiuj). We will focus on the mexican hat wavelets given by: 

g-lx-?lV2i?^ (15) 

From the wavelet transform of a function f{x) one can obtain the wavelet coefficients : 

w{b,R) = J dx^{x;b,R)f{x) (16) 

and if f{x) is Gaussian, w{b,R) will be Gaussian as well. We will here use the CMB 
temperature fluctuation field as the f{x) function. We will implement the non-Gaussianity 
test in wavelet space as we did for the Minkowski functional: 

• we generate a set of Gaussian CMB maps for calibration and a set of non-Gaussian 
maps for testing 

• we cut tangent planes (12° x 12°) 

• we calculate for each plane the coefficients w{b, R) 

• we evaluate the skewness of the wavelet coefficients for each sky. 

• finally, using the skewness of the wavelets from the Gaussian maps, we define the one 
and two a detection levels as described above for the other methods. 



^(x; b, R) 



(27r)i/2i? 



\x 



R 



VII. COMPARISON AND COMBINED TEST 



This section aims at comparing the different methods described above. Applying the 
methods on the same maps , we will first compare the number of detections. We will in this 
paper use the non-standard infiationary model described in Q| which has the non-linear 
coupling parameter /at^ as a measure of the strength of non-Gaussianity. We generated 100 
'observed' skies with Planck-like noise (LFI 100 GHz), beam 20', pixel-size ~ 6' (Nside 512 
in Healpix language), using a pure Sachs- Wolfe spectrum with /nl values of 300 and 100. 
In figures |21 ID and El one can see the behavior of the Minkowski functionals in the presence 
of a non-zero J'nl. In tables H] and |n] we show the rejection rates. As the first Minkowski 
functional was giving the best results in this tests, we will focus only on Mq for this kind 
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of non-Gaussianity. In table IIIII we list the results of the wavelet test on the same maps 
(for the wavelets we used the parameter R = 22.5'). For individual results of the empirical 
process test, we refer to 3^ . 



TABLE I: The quantile levels which determine the one and two a detections and the rejection rates 
for non-standard inflationary models with /at/, = 300 



quantile detection limits 


lo 


la 


3.6 10-3 


h 


la 


0.286 


h 


la 


33.2 




2a 


6.3 10-3 




2a 


0.524 




2a 


49.7 


Rejection rates % 




la 


100% 


h 


la 


100% 


h 


la 


100% 




2a 


100% 




2a 


71% 




2a 


84% 



TABLE II: The quantile levels which determine the one and two a detections and the rejection 
rates for non-standard inflationary models with fNL = 100 



quantile detection limits 


lo 


la 


3.6 10-3 


h 


la 


0.286 


h 


la 


33.2 




2a 


6.3 10-3 




2a 


0.524 




2a 


49.7 


Rejection rates % 


lo 


la 


68% 


h 


la 


52% 


h 


la 


57% 




2a 


35% 




2a 


8% 




2a 


8% 



Table IIVI shows the number of detections at the different levels, using Mq and the 
empirical process method on maps with f^i = 100 and /nl = 300. The power of the two 
procedures appears very close. However, analyzing the individual maps, we find that only 
one third of the maps detected at 2a are the same for the two tests. This leads to the idea 
of implementing a combined test. 
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TABLE III: wavelets test: the rejection rates at one and two a for non-standard inflationary models 
with /iVL = 100 



fNL = 100 


Confldence Level 


Rejection Rate 


1 a 

2 a 


89% 
57% 



TABLE IV: COMPARISON AND COMBINED TEST 



TEST 


emp. proc 




Mo 




comb 


Inl 


100 


300 


100 


300 


100 


1 a 


50% 


95% 


68% 


100% 


74% 


2 a 


29% 


87% 


35% 


100% 


35% 



For the combined test we suggest to use an indicator consisting of Jq from the Minkowski 
functionals and fcmax from the empirical process. We chose to normalize the Jq and fcmax so 
that they both have mean zero and variance one, using Monte-Carlo simulations of Gaussian 
maps. In this way, the two values can be averaged: 

X = wJo + W2fcmax, (17) 

where 

io = (/o - (/o))/v^a^r^ao)^, (18) 

and 

kmax (^max (^max))/ \/(^max) (^max) • (l"^) 

Here () means mean value taken over 100 Gaussian simulations. The weights wi and W2 
were chosen proportionally to the power of each procedure. Using the rejection rates for 
Inl = 100 in table HIl we arrived at Wi ~ 0.6 and W2 ~ 0.4. Of course, the threshold 
value for x needs to be evaluated anew. In figure (jHl) we plot the distribution of x for the 
Gaussian and f^L = 100 non-Gaussian maps. 
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TABLE V: The results of the combined test using all three methods: Empirical process, Minkowski 
functionals and wavelets for /nl = 100 



fNL = 100 


Confidence Level 


Rejection Rate 


1 a 

2 a 


100% 
78% 



By a similar motivation, it is natural to combine also the wavelet method into a single 
procedure; it is easy to see (tables ITTl and that the detection rate at 2a is about two times 
higher for the wavelets, warranting the wavelet coefficient a very high weight in the combined 
analysis. Inspecting table lIVl we detect a moderate improvement in the detection rate when 
combining Minkowski functionals and empirical process. The result of the empirical process 
+ Minkowski functionals + wavelets combined test is shown in table |3 A significant 
improvement of the number of detections at both confidence levels is evident; note that 
the combined procedure is to some extent model dependent, as the weights we used were 
tabulated from specific non-Gaussian models. 



VIII. POINT SOURCES 



We also simulated maps with point sources (with noise and beam as given above) to 
compare the power of the methods on a different kind of non-Gaussianity. We generated a 
toy model of point sources with a distribution building on formula (1) in 4l| and formulae 
(1) and (2) in In table IVTl we show the results for the Minkowski functionals. We see 
immediately that for this kind of non-Gaussianity, the first Minkowski functional is not 
sensitive, whereas the other two functionals show a good rejection rate. This suggests that 
we might be able to discriminate between these two types of non-Gaussianity. The first 
Minkowski functional can be used to trace primordial non-Gaussianity with little infiuence 
from the point sources. On the other hand, the presence of point sources will show up in the 
second and third Minkowski functionals, which are only weakly infiuenced by primordial 
non-Gaussianity. In Figure (jZj) we show the shape of the Minkowski functionals in the pres- 
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FIG. 6: Histogram of the combined estimator x (empirical process and Minkowksi functionals) 
for Gaussian realizations (solid line) and non-Gaussian realizations with /nl = 100 (dotted line). 
All realizations have Planck-like noise and beam 20'. The shaded areas represent the la and 2a 
detection limits. 

ence of point sources. Note that the deviations from the Gaussian mean are different than 
in the case of primordial non-Gaussianity (Figure 0}. The point sources manifest themselves 
mainly as an offset in Mi and M2, consistent with what was observed for weak lensing 
j^]. However as seen in figure (0) for non-standard inflation the curve has a particular shape. 

One could imagine combining Mi and M2 in order to strengthen the power of the test, 
similarly to what we have done above. However, it turns out that the maps detected by M2 
are contained within the maps detected by Mi, so that there is no additional information 
in combining the two estimators. 

For the empirical process method, there was no detection for any of the tests, univariate, 
bivariate or trivariate. The numbers kmax obtained for the maps with point sources were 
consistent with those for Gaussian maps. Also for the wavelet test, the number of detections 
was very small. Note that wavelets can be useful for detection of bright point sources j^l ) 
but in our source model these were excluded. This might suggest that the empirical process 
methods and the skewness of the wavelets can be used to probe primordial non-Gaussianity 



16 



Gajssian vs Non-Gaussian, (point sources) 




Gaussian vs Non-Gaussian, (point sources) 



-3 -2 -1 



Gaussian vs Non-Gaussian, I^point sources) 




threshold (i') 



FIG. 7: Minkowski functionals (divided by the mean) averaged over 200 realization of Gaussian 
maps (solid line) and 100 non-Gaussian maps (squares) with point sources, Planck noise and beam 
20'. The shaded bands show the la deviations of the Gaussian realizations and the error bars show 
the la deviations for the non-Gaussian maps. 



without confusion from point sources, making the combined test presented above more 
robust. 
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TABLE VI: The quantile levels which determine the one and two a detections and the rejection 
rates for maps contaminated by point sources 



quantile detection limits 


h 


la 


1.1 10-3 


h 


la 


0.415 


h 


la 


207.6 




2a 


1.7 10-3 




2a 


0.770 




2a 


328.3 


Rejection rates % 


h 


la 


36% 


h 


la 


81% 


h 


la 


75% 




2a 


8% 




2a 


45% 




2a 


33% 



IX. COMMENTS AND CONCLUSIONS 



We have compared three methods for detecting non-Gaussianity in observations of the 
cosmic microwave background. The methods were applied to non-Gaussian maps with two 
different kinds of non-Gaussianity, primordial non-Gaussianity with a varying /^vl and 
point sources. It is important to note that the non-Gaussian maps used in this article 
were generated taking into account only the Sachs- Wolfe effect. In future work we will 
study maps where the full radiative transfer equations have been applied. For the time 
being, we stress that our results are broadly consistent with the power of the procedures 
adopted for WMAP data analysis. More precisely, a detailed comparison is unfeasible, as 
we are assuming a simplified non-Gaussian model (no radiative transfer) and Planck LFI 
like noise and beam. Moreover we are restricting the analysis to the first 500 multipoles 
. Broadly speaking, however, Monte-Carlo simulations suggest that a value of /atl about 
100 represents the lower limit that can be detected at a 2a level, by usin g a combined 
procedure: this seems consistent with the value fjqi = 139 reported in ([ll|). As an 
estimator of non-Gaussianity for the Minkowski functionals we have introduced a statistic 
which gathers the information from different thresholds feqllO|). The estimator for the 
empirical process method is kmax, the maximum value of the function K{a,r) obtained 
from a given map (eqUTj) using the trivariate test. For the test in wavelet space, we use the 
skewness of wavelet coefficients. 
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For the primordial non-Gaussianity, Mq and the empirical process method have a similar 
rejection rate, whereas Mi and M2 showed less power. On the other hand, for the maps 
with point sources. Mi and M2 gave the best results, whereas Mq and the empirical process 
method had no rejections. The fact that the first Minkowski functional shows little power 
in the presence of point sources is hardly surprising. Indeed this statistic depends only on 
the pixel by pixel temperature values and hence is not at all affected by discontinuities in 
the map. The converse is clearly true for the other Minkowski functionals, which are sen- 
sitive to the local morphology of the maps. For the empirical process, we simply note that 
spikes in real space are erased in harmonic space. It is also important to stress that our 
results depend heavily upon the nature of non-Gaussianity; in particular some preliminary 
exploration of Monte-Carlo evidence from non-physical toy models suggests that the power of 
these procedures need not be close, in general. This strengthens the case for (weighted) mul- 
tiple/combined procedures; the combined procedures seem to show a marked improvement 
in the power of the test. However, the pixel-, harmonic- and wavelet-space methods, despite 
carrying complementary statistical information, should not be viewed neither as orthogonal 
nor as independent, so that some care is needed when merging them into a single statistic. 
In any case, the fact that different methods detect different kinds of non-Gaussianity can be 
viewed as an advantage, in the sense that, for instance, primordial non-Gaussianity can be 
detected without confusion from point sources. 
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